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We propose a new model suitable for a nonequilibrium molecular dynamics (MD) simula- 
tion of electrical conductors. The model consists of classical electrons and atoms. The atoms 
compose a lattice vibration system. The electrons are scattered by electron-electron and 
electron- atom interactions. Since the scattering cross section is physically more important 
than the functional form of a scattering potential, we propose to devise the electron-atom 
interaction potential in such a way that its scattering cross section agrees with that of 
quantum-mechanical one. To demonstrate advantages of the proposed model, we perform 
a nonequilibrium MD simulation assuming a doped semiconductor at room or higher tem- 
perature. In the linear response regime, we confirm Ohm's law, the dispersion relations and 
the fluctuation-dissipation relation. Furthermore, we obtain reasonable dependence of the 
electrical conductivity on temperature, despite the fact that our model is a classical model. 
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1. Introduction 

In statistical mechanics, construction of nonequilibrium statistical mechanics has been 
attempted for long years. Whereas in the linear nonequilibrium regime the linear response 
theory was established in the 1950s, 1 ) in the nonlinear nonequilibrium regime the properties 
of nonequilibrium states far from equilibrium are still poorly understood. 

When trying to investigate such nonequilibrium states with analytical approaches, one 
runs into difficulties of solving the equations of motion analytically. A promising approach is 
the nonequilibrium molecular dynamics (MD) simulation, which came into sight in 1950s. 2 ) 
In this approach, assuming an appropriate microscopic model, one numerically solves the 
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equations of motion of all constituents, obtains the values of macroscopic variables of interest, 
and thereby understands properties of nonequilibrium states. 

Transport phenomena are the most important in nonequilibrium statistical mechanics. 
For heat conduction, many studies by the nonequilibrium MD simulations were reported. For 
example, the Fourier law was successfully obtained. 3 ) For charge or mass transport, in contrast, 
the nonequilibrium MD simulation was successfully performed only recently by Yuge, Ito and 
Shimizu in 2005. 4 ) Their model, which we call the YIS model, consists of three types of classical 
particles, which imitate electrons, phonons, and impurities. These particles interact with each 
other via short-range interaction potentials. Using the YIS model, it was subsequently shown 
for states far from equilibrium that response becomes strongly nonlinear, 4 ) the long-time 
tail is significantly modified, 5 ' 6 ) the fluctuation-dissipation relation (FDR) is significantly 
violated and universal excess noise appears. 7 ) It was further shown that the sum rules and 
the asymptotic behaviors, which were recently derived in refs. 8-10, of response functions of 
NESSs are indeed satisfied. 8 ) 

Although the YIS model have enjoyed such great successes in treating fundamental prop- 
erties of NESSs, it has two problems to treat more general properties. First, we have found 
that the YIS model is not suitable for treating macroscopically inhomogeneous conductors, n ) 
because in NESSs of such conductors the phonons are pushed by the electrons away from the 
conductors, which is physically unrealistic. 12 > Second, the YIS model is inconvenient for the 
analysis of the temperature dependence of the electrical conductivity, because the number of 
phonons has to be changed as a function of temperature. These difficulties arise because the 
phonons are treated as classical particles. 

In this paper, we propose a new model which resolves these difficulties. In this model, we 
represent a phonon system by a classical lattice vibration system. Each electron is scattered 
by lattice vibrations and by other electrons (and possibly by impurities). As an illustration of 
advantages of the present model, we will show by MD simulations that not only nonequilibrium 
properties at each temperature but also the temperature dependence of the conductivity can 
be analyzed in a natural way. 

This paper is organized as follows. We describe the essential elements to treat electrical 
conduction in Sec. 2. The new model, which includes all the essential elements, is presented in 
Sec. 3. To illustrate advantages of the model, we perform an MD simulation in Sec. 4 assuming 
a doped semiconductor at room or higher temperature. The results will be presented and 
discussed in Sec. 5. We devote the last section to the summary. 
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2. Essential Elements to Treat Electrical Conduction 

Before defining the model, we here present the essential elements to treat electrical con- 
duction. 

A typical experimental situation is as follows. An electrical conductor (e.g., a doped semi- 
conductor) is put on a sample holder whose temperature is kept constant. The both ends of 
the conductor are connected to a battery through conducting wires. By inserting an ampere 
meter to the wire, one measures the conductance by the two- or four-terminal method. 

From a microscopic viewpoint, electrical conduction is explained as follows. Energy is 
supplied to the conductor from the battery, and global motion of the electrons is induced, 
which results in a finite electrical current. The supplied energy is transferred to the lattice 
through the electron-lattice interaction, and dissipates as the Joule heat into the sample holder 
through the lattice-bath coupling. A balanced state in which the supplied energy equals the 
Joule heat is a nonequilibrium steady state (NESS). 

Therefore, a model for electrical conduction should have the following elements: (a) an 
interacting many-electron system, (b) an energy dissipating system, such as a lattice vibration 
system and a bath for the lattice, (c) objects violating the microscopic translational invariance 
to define the rest frame, such as impurities, fixed walls, and the lattice of atoms (d) interactions 
among these constituents, and (e) an external force to drive electrons. 

3. Proposal of Model 

We propose a model that includes all the elements in the previous section. From now on, 
we explain the following components one by one: a many-electron system, a lattice vibration 
system, a proper interaction between these two systems, impurities, an external electric field 
and the boundary conditions. 

3.1 Many- Electron System 

We consider the regime where the electron temperature k-gT c is higher than the Fermi 
level £f (the chemical potential at zero temperature): k-QT e > eF- In this regime, we can treat 
the electrons as classical particles. 

The kinetic energy of the classical many-electron system is given by 



where m is the effective mass of an electron, and V{ is the velocity of the electron labeled by 



We assume a conductor which is translation-invariant macroscopically. Then, the electron- 




(1) 



i (i = l,2,...,iV e ). 
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electron interaction potential U ee is reduced to a short-range potential, because long-range 
effects of the Coulomb potential are screened. It is expected that the detailed form of short- 
range potential U ee has no significant effect on transport properties. Hence, among many 
possible forms of short-range potentials, we here take the following simple one: 

u ™ = Y, Y, ( l - |r ^~ rjl ) Q (^c - In - rj\) , (2) 

i j>i ^ ee J 

where Q ee is a positive constant, and r,- L is the position of the zth electron, and / ee is the 
interaction range, and 6 is the step function. The factor @(/ C c — \fi — ensures that 
pairs of electrons interact only when the electrons are within the distance l ee . The factor 
(1 — \ri — rj\/l cc ) 4 is introduced to weaken the singularity of U cc at |rj — Vj\ = l cc . When the 
constant Q ee is much larger than ksT c , U cc can be regarded as a hard core potential. 



3.2 Lattice Vibration System 

We consider the case where the temperature is higher than the Debye frequency: k^T > 
huii). In this case we can treat the quantum lattice vibration system as a classical lattice 
vibration system. It consists of many classical atoms, whose kinetic energy is given by 

a 

where M is the mass of an atom, and V a is the velocity of the atom labeled by a (a = 
l,2,...,JV a ). 

Suppose that the atoms compose a square lattice, whose lattice constant at mechanical 
equilibrium is denoted by I. We assume that each atom is connected to its nearest-neighbor 
atoms by nonlinear springs. The atom-atom interaction potential C/ aa is taken as 




where the summation is taken over all pairs of nearest neighbors. R a is the position of the 
crth atom, lo is the natural length of each spring, and K and G are positive constants. Since 
the nonlinear terms are included in U aSL , the lattice system is expected to be chaotic if iV a 
is large enough. It is therefore expected to have good statistical properties. Furthermore we 
impose the condition 

I > lo , (5) 
to prevent the zero-frequency angular modes. 
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We also assume the following self potential: 




a 



(6) 



where K' and G' are positive constants, and R a o represents the mechanical equilibrium po- 
sition of the ath atom. This potential stabilizes the positions of the atoms around their 
mechanical equilibrium positions more rigidly. 

We impose the condition that the lattice doesn't melt; 



where 5 is the standard deviation of the position of an atom: 5 = a/ (\R a — Rao\ 2 ) ({A) 



Hence, even if momentum of an electron (~ ^/k-BT c m) is fully transferred to an atom, its effect 
on the vibrating motion of the atom is quite small. 

3.3 Interaction between Many-Electron System and Lattice Vibration System 

In this subsection, we propose a classical-mechanical form of the electron-atom interaction 
potential ?7 ca , which reproduces a quantum-mechanical scattering cross section. 

3.3.1 Electron-phonon scattering in a quantum-mechanical system 

Before discussing a classical form, we briefly review the theory of electron-phonon scat- 
tering in solids. 13 ^ 14 - > 

Consider the single-electron potential produced by the lattice. When the lattice is perfectly 
periodic, the potential (denoted by Uo(r)) is also periodic: Uo(r) = Uo(r + n), where n is 
the lattice vector. In this case we can incorporate the effect of the lattice potential on an 
electron by replacing the electron momentum with the crystal momentum, and by modifying 
the electron's dispersion relation. When the electron energy is much smaller than the band 
width, we can employ the effective-mass approximation. 

When the positions of the atoms are displaced from the periodic structure, the lattice 
potential (denoted by Ud(f)) is non-periodic, and the difference, V = Ud{r) — Uo{r), causes 
the electron-phonon scattering. 

The scattering cross section s quajltam in the Born approximation is proportional to 
|(/|y|i)| 2 where \i) and |/) are the states before and after the scattering, respectively. Its 



I > 5 , 



(7) 



denotes the averaged value of A), and 5 ~ ^Jk^T jK when G = K' = G' = 0. 
We note that the mass of an electron is much smaller than that of an atom; 



m < M . 



(8) 
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magnitude depends on the phonon mode. 

For the scattering by the acoustic modes in the long wave length limit, for example, 
V ~ E{\7 ■ d(r), where E% is a constant, and d is the displacement of an atom from its position 
on the perfect lattice. In the Fourier space, Vq oc q ■ d{q). Therefore s qua ntum °c |(/|^|i)| 2 oc 
|V • d(r)\ 2 oc |<i||(q)| 2 , where d\\(q) denotes the component of d(q) parallel to q. 

3.3.2 Electron- atom scattering potential in the classical-mechanical model 

Since we are considering the case where k&T c > ep, quantum interference effects are weak, 
and the electrons may be treated as classical particles. Let us consider the scattering potential 
U ea for the classical electrons. 

Generally, to model a quantum-mechanical system using a classical-mechanical system, 
the scattering cross section is much more important than the functional form of a scattering 
potential. 15 ^ 

In the present case, if one employed that the classical potential U ea that has exactly the 
same form as the quantum potential U eSL , then the classical electrons would be scattered by U ca 
even when the atoms are not displaced from the lattice points. This effect sharply contradicts 
the nature of the quantum system, in which electrons are not scattered (except for Umklapp 
processes) by a perfect lattice of the atoms. 

We therefore propose that U ea should be devised in such a way that the classical scattering 
cross section ^classical imitates the quantum scattering cross section s qU antum well. Specifically 
when we consider the case where s qU antum oc \d\ 2 , U ca should lead to ^classical oc 5 2 for each 
atom. If U cai satisfies this condition, it is expected that the detailed form of U ca is irrelevant 
to transport properties. 

As an example of such U ca in a two dimensional system, we here take 16 ) 



X 6 7e 




(9) 



I 

where Q ea and 7 ea are positive constants which characterize the magnitude and the interaction 
range of U ca , respectively. Due to the step function 0, the interaction range (for an electron) 
of this potential is ~ 7 ca |-Rg — Rpo\ 2 /l (I is the lattice constant), which is proportional to 
the square of the displacement \Rp — Rpo\ of the atom. Since s c i ass ; ca i is proportional to the 
interaction range in a two dimensional system, we have s c i aS sicai oc 5 2 , as required. The factor 
(l — l\r a — Rp\/"fea.\Rp — Rpol 2 ) 4 is introduced to weaken the singularity at |r Q — Rp\ = 
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7 ea |i?/3 — Rpo\ 2 /l. Although this factor is singular with respect to Rp at Rp = Rpo, this 
singularity is removed by the extra multiplicative factor \Rp — Rp /l\ 8 . 

For a three dimensional system, ^classical is proportional to the square of the interaction 
range. We can therefore take U ea , for example, as 



u °» = EE^ c 



R/3 — R/30 



X 6 7l 



I 

Rp — Rpo 



l 



Re 



I 



7ea|-R/3 - Rpo\ 

t - Rf3 



I 



(10) 



This gives s c iassicai oc 5 2 , as required, in a three dimensional system. 



3.4 Impurities 

To violate the microscopic translational invariance and thereby define the rest frame of 
equilibrium states, we add a random potential. We represent the random potential as the sum 
of short-range potentials produced by impurities. The impurities are modeled by particles 
which are fixed at random positions and interact with the electrons and the atoms through 
the short-range potentials. Among many possible forms of the short-range potentials, we 
here take the following forms for the electron-impurity potential U c i and the atom-impurity 
potential t7 ai : 

^ei = E E (l " ^T-^) 4 " - , ("J 

i k * el 

U ai = E E 3- ( 1 " lRa ~^ ) A Q(ki ~ \Ra ~ &|) , (12) 
a k V £ ai / 

where £ k is the position of the kth impurity (k = 1, 2, N{). Q e i and Q a i are the positive 
constants. l c \ and l & \ are the interaction ranges. 

3.5 Boundary Conditions 

To sum up, we have obtained the following Hamiltonian: 

H = K c + K a + U cc + U ca + c7 aa + U & + U e \ + C/ai . (13) 

In addition, an electric field in the x-direction E is applied to all electrons, and each electron 
experiences a force F = eE. 

Hereafter, we consider two-dimensional systems, employing eq. (9) in §3.3.2 as the electron- 
atom scattering potential. A schematic picture of the model is shown in Fig. 1. 

The boundary conditions are imposed as follows. The boundaries in the x-direction are 
the periodic boundaries for all particles. The boundaries in the y-direction for electrons are 
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potential walls, which simulate the walls at the boundaries of the conductor. Those for atoms 
are thermal walls, which are located away from the potential walls for electrons by a dis- 
tance Lft (see Fig. 1). The thermal walls simulate the thermal contact with a sample holder 
whose temperature is kept at T. Through these thermal walls, heat is transferred outside the 
conductor, and a NESS can be realized in the conductor. 

Note that the present model is a pure mechanical model, except for the thermal walls for 
atoms. 




E 

Fig. 1. A schematic picture of the two-dimensional model of a macroscopically homogeneous con- 
ductor. The large black circles and the small black circles represent electrons and impurities, 
respectively. The gray circles represent atoms, which compose a square lattice. By an electric field 
E, the electrons move from the left to the right in average, interacting with the atoms and the 
impurities. 



4. Application to Doped Semiconductors 

To illustrate advantages of the proposed model, we apply the model to doped semicon- 
ductors, and perform an MD simulation. 
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4-1 Doped Semiconductors 

The proposed model can be applied to doped semiconductors at room or higher tempera- 
ture (but lower than the melting temperature). We explicitly show this for the case where the 
electrons are confined in a two-dimensional plane (e.g., in a quantum well or inversion layer) 
by showing that the conditions of Sec. 3 are satisfied. 

Since the electron temperature T e in NESSs is higher than the bath temperature T, k^T c > 
k^T > fcB?r 00m ~ 26meV. On the other hand, ep — irh 2 n e /m where n e is the electron density. 
For n c — 5.0 x 10 15 m~ 2 , for example, eF — 6.5meV in doped silicon (Si) semiconductors, 
whereas eF — 18meV in doped gallium arsenide (GaAs) semiconductors. In both cases, the 
condition fee^e > €f is satisfied. 

The condition k-^T > Hujd is approximately satisfied, because in Si and GaAs crystals, the 
Debye temperature is 645K and 633K, respectively; hence k-B,T ~ Swd. 

The condition I S> S ~ ^/ksT/K is also satisfied, because we are considering the case 
where the temperature is less than the melting temperature. (In Si and in GaAS crystals, the 
melting temperature is 1687K and 1511K, respectively.) 

Generally, at room or higher temperature, electron-atom scattering occurs more frequently 
than electron-impurity scattering. Hence, we will take the model parameters in such a way 
that 

Tea < T ei (14) 

is satisfied, where r ea is the mean free time between electron-atom collisions, and r c ; is that 
between electron-impurity collisions. 

Therefore, our model is applicable to typical doped semiconductors. 

4.2 MD Simulation 

We analyze electrical conduction in doped semiconductors by an MD simulation of the 
proposed model. We set the lattice constant Z, the effective mass of an electron to, the electric 
charge e, and a certain reference energy to unity. 

As illustrated in Fig. 1, electrons and impurities are confined in —L x /2 < x < L x /2 and 
—Ly/2 < y < Ly/2, whereas atoms are confined in —L x /2 < X < L x /2 and —(L y /2 + Lb) < 
Y < Ly/2 + L b . 

At an initial time (t = 0), we put atoms at their mechanical equilibrium positions, and 
arrange electrons and impurities at random positions so as not to contact with each other. 
The initial velocities of electrons and atoms are given by Maxwell distribution of the bath 
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temperature T. In other words, the initial state is closed to an equilibrium state. 

We set the numbers of electrons (N e ), atoms (iV a ), and impurities (iVj) by considering the 
density of a real doped semiconductor. In our simulation, we take the density of electrons high 
(N e /N a ~ 1/3) to shorten the computational time, although in a usual doped semiconductor 
the density of electrons are rather low (N e /N a <C 1). 

An external electric field E is applied for t > 0. Electrons and atoms move according to 
Newton's equations of motion. We use Gear's fifth-order predictor-corrector method to solve 
the equations of motion. We also use the Linked Cell method 19 ) in calculating the force due 
to U cc , U ea , U ci and C/ ai . 



We set k-gT = 1, except when we investigate the temperature dependence. We set the 
parameters as follows: Q cc = 10000, l ee = 1/2, M = 10, Q ca = 10 10 , 7 ca = 1, Qci = 10000, 
Z e i = 1/12, Qai = 10000, Z ai = 1/4, K = 100, G = 800, K' = 100, G' = 800 and l = 3/4. Here, 



electrons' flow too much. For these values of the parameters, effects of the nonlinear terms 
are small, because K5 2 /2 > G5 4 /A and K'5 2 /2 > G'5 4 /A. These parameters satisfy the 
conditions in the previous sections: I > lo, I S> S, M S> m, and r ca <C r e ; (for the particle 
densities assumed in the following MD simulations). 

In the x-direction, we take the periodic boundary condition. In the y-direction for electrons, 
the boundaries are the potential walls. The wall potential is taken as U = Q ee (w/l ee ) 4 , where 
w is the penetration depth of an electron into the wall. For atoms, the end atoms around 
\y\ = Ly/2 + Lf, are connected to hard walls which are located at a distance I away from the 
mechanical equilibrium positions of the end atoms (i.e., at \y\ = L y /2 + Lb + 1). When an end 
atom comes in contact with a thermal wall (located at \y\ = L y /2 + Lb), 20 ) the atom loses 
its memory and is reflected back with a random velocity 21 - 1 whose probability distribution 
function is given by 



macroscopic behaviors of the system are almost independent of the time-step width if the 
width is taken smaller than this value. 

We calculate the electrical current as follows. The x-component of the electrical current 
density at time t is given by j(x,y,t) = J2i e v f(t) ~ x i)Hv ~ Vi)i where vf(t) denotes 
the x-component of the velocity of the zth electron. The current in the x direction at x 
is given by I(x,t) = f dy j(x,y,t) = J2i e vf(t)S(x — xi). When measuring the current in 



we have taken Q eil (5 /I) 8 > k^Te ~ ksT to estimate r ea easily, and l e i small not to impede the 




(15) 



In the main simulations, we fix the time-step width to be 10 4 . We have confirmed that 
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an experiment, one usually measures the current averaged over some finite length L in the 
x-direction: I(t) = f dx I(x,t)/L = e ^ vf (t)/L . In this simulation, we take L = L x , i.e., 

I(t) = e £ (16) 

ig conductor 

5. Results of MD simulation 

5. 1 Snapshot 

Figure 2 is a snapshot of a NESS at E = 0.005 for a small-size system with L x = 30, L y = 
10, Lb = 1,N C = 150, iV a = 390, iV; = 10. We observe that atoms (gray circles) compose a 
well-defined lattice, whereas electrons (large black circles) distribute rather randomly. Similar 
snapshots (not shown) are obtained for larger-size systems used in the following simulations. 



Ly/2 + Lb 



« © »■ • •*•••**• Ly/2 

. . . *? v. .......... * f» 



r . # . • • . . 
*• .. j . . . 

• • • • .• • • • f^< 

# .•.»». . . « . 



• • • 
, * ••• • • 



••• t — ** • 

^ • • • • • 

a» • * • • . *• •'• • « # « • • ••••••• 

•:. •••• :•* , ■ • : • * 
. . * • v • •• • •• •• • ••••• 

• • •••• #■ 





Lx/2 







Lx/2 
\ 





Lj/2 



E 



Fig. 2. A snapshot of a NESS in a conductor at E = 0.005. The large black circles, the gray circles, and 
the small black circles represent the typical interaction ranges of electrons, atoms, and impurities, 
respectively. 



5.2 Realization of a NESS 

We first investigate whether the system reaches a NESS in the presence of an external 
electric field. A NESS is said to be realized if every macroscopic variable fluctuates around 
a constant value and if the fluctuation is relatively small for a large-size system. 8 ' 9 ) We here 
examine two macroscopic variables: the electrical current in the x-direction, I(t), and the total 
kinetic energy, K(t). To smear out high-frequency components, which will make the figures 
too busy, we plot I(t) and K(t) after averaging over time intervals At = 6. 

Figure 3 depicts the time evolution of I(t) for E = 0.005 for a large-size system with 
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L x = 100, L y = 36, L b = 2, N c = 1500, iV a = 4100, N { = 100. For t > 800, we see that I(t) 
fluctuates around the long-time-averaged value (represented by the solid line in the figure). 
Although the fluctuation might look rather large, its magnitude is reasonable for this size of 
system. In fact, the fluctuation is of the same order of magnitude as that in an equilibrium 
state (E = 0), and the latter is consistent with the fluctuation-dissipation relation as will be 
shown in §5.5. 

Figure 4 depicts the time evolution of I(t) for larger E (E = 0.04). Comparing Figs. 3 
and 4, we see that the relative fluctuation for large E is smaller than that for small E. 




time t 

Fig. 3. Time evolution of the electrical current I(t) at E = 0.005. The solid line shows the time- 
averaged value of I(t) during 1000 < t < 2000. 

Figure 5 depicts the time evolution of K(t) = K c (t) + K a (t) for E = 0.005, where K e and 
.Ka are the kinetic energies of electrons and atoms, respectively. For t > 1000, K(t) fluctuates 
around the long-time-averaged value, and the fluctuation is small enough. Furthermore, K(t) 
is fitted well to the function, 

K(t) = K(0) + a(l - e"* /r - lax ), (17) 

which is represented by the solid line in the figure, where K(0) = 3800, a = 2700, r rc i ax = 280. 
This indicates that influence of the initial state becomes very small for t > 1000 (because 



12/21 



J. Phys. Soc. Jpn. 



Full Paper 




CD 

' 1 1 1 

500 1000 1500 2000 

time t 

Fig. 4. Time evolution of the electrical current I(t) at E = 0.04. The solid line shows the time- 
averaged value of I(t) during 1000 < t < 2000. 

exp(- 1000/280) ~ 0.028). 

From these observations, we conclude that the system reaches a NESS in the presence of 
an external electric field after a sufficiently long time. Hence, in the analysis below, we will 
use the time-averaged values calculated for t > 1000 for the average values of macroscopic 
variables in a NESS. 

5.3 Dependence on Electric Field 

Next, we investigate the dependence of the electrical current I on the electric field E. 
We calculate long-time averaged values of I independently for five systems with different 
impurity configurations, and calculate the average and standard deviation of the long-time 
averaged values. This is because we want to investigate the conductivity a = I/L y E in a 
sufficiently large-size system. In an MD simulation, we can calculate the conductivity only in 
an insufficiently large-size system. We therefore take five impurity configurations, and estimate 
the conductivity in a sufficiently large-size system. 

The result is shown in Fig. 6, where the long time averaged value of I(t), denoted by /, 
is plotted against E. We see that conductivity is nearly constant for E < 0.003. Therefore, 
Ohm's law I <x E holds for small E (E < 0.003), whereas the nonlinear response is observed 
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Fig. 5. Time evolution of the kinetic energy K(t) as a function of time at E = 0.005. The solid line 
is the fitting curve whose functional form is given by cq. (17). 

for larger E. The standard deviations are very small, which implies that the dependence on 
the impurity configuration is relatively small. 

We also investigate the dependence of the electron temperature fcsT,? = m(vf — (vf) 2 ) on 
the electric field E. As plotted in Fig. 7, we find that 



holds for E < 0.005 where 6 is a positive constant (whereas k B T* increases more slowly for 
E > 0.005). This curve reflects the fact that the Joule heat is proportional to E 2 in the linear 
response regime. To see this point, in Fig. 8, we plot the dependence of k B T£(E) on the Joule 
heat E ■ I. For small E, it is clearly seen that k B (T^(E) — T£ (0)) oc E I, as expected. For 
larger E, however, Tg increases more slowly. This means that energy transfer from electrons 
to phonons becomes very fast for larger E. 

5.4 Dispersion Relations 

We study whether the dispersion relations (Kramers-Kronig relations) of the complex 
conductivity ct(oj) holds in this model. To obtain a(oj), we apply an AC electric field E(t) = 
Eq cos(otf), measure the current during 1000 < t < 10000, and calculate the real and imaginary 
parts of the conductivity a{oS) = I{oS) / 'L y E{uj), where ~ denotes the Fourier transform. We 



k B T* ~ k B T + bE 2 



(18) 
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Fig. 6. Average values of the electrical current I plotted against the electric held E. The error bars 
indicate the standard deviation over five impurity configurations. The dotted curve represents the 
linear response. 
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Fig. 7. Average values of the electron temperature k^T^ against the electric field E. The error bars 
indicate the standard deviation over five impurity configurations. The dotted curve represents eq. 
(18) with b = 2.2 x 10 4 . 
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Fig. 8. Average values of the electron temperature k&T* against the Joule heat E ■ I. The error bars 
indicate the standard deviation over five impurity configurations. The dotted line represents the 
linear relation: k B (T*(E) - T e x (0)) oc E ■ 7. 



set Eq = 0.001, which corresponds to the linear response regime. 

For a single impurity configuration, we calculate <j(oj) at each co by averaging the results 
for 10 samples, which have different initial distributions of electrons and atoms, as well as 
different random numbers for the thermal walls. 

The result is plotted in Fig. 9. To check the validity of the dispersion relations, 
duJ 1 V , , „ , , , f 00 duo' V 



Re 



-r Ima(u') , Ima(oj) = - / — — Rea(oj') , (19) 

.oo IT OJ' -UJ J_ 00 TT CO 0J 



we need to perform the Hilbert transformation. However, the transformation is generally 

difficult because it requires knowledge of Re ct(uj) and Im a (to) over a very wide range of u. 

Instead, we fit the data to the Drude formula: 

u i \ n ee 2 r n e e 2 ujT 2 

m(l + (ojtY) m(l + [UJT) Z ) 

which satisfy the dispersion relations. The fitting parameter is r. We see that the data in 

Fig. 9 are well fitted by eq. (20) with r = 59. Therefore, we conclude that the dispersion 

relation holds. 
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Fig. 9. A semi-logarithmic plot of the real and imaginary parts of the complex conductivity <r(w) 
against the angular frequency u>. The error bars indicate the standard deviations over 10 samples. 
The solid and dotted lines are the real and imaginary parts of the fitting curves, eq. (20). 



5.5 Fluctuation- Dissipation Relation 

We also examine whether the fluctuation-dissipation relation (FDR) holds in this model. 
Since the complex admittance is given by a(uj)L y /L x , the FDR is expressed as 

= 2k B T Re o-(u>) ■ ^ , (21) 

where Si(u) is the spectral intensity of the electrical current in the equilibrium state (E = 0), 
and a(u)) is the complex conductivity in the linear response regime (Eq = 0.001). We measure 
the electrical current during 1000 < t < 10000. For a single impurity configurations, we obtain 
the values of the left-hand side (LHS) from 20 samples, and the values of the right-hand side 
(RHS) from 10 samples. 

The results are plotted in Fig. 10. We see that the LHS and the RHS agree well at each 
value of uj. Therefore, we conclude that the FDR holds in this model. 

5.6 Temperature Dependence of Electrical Conductivity 

We investigate the temperature dependence of the DC electrical conductivity a(T) in the 
linear response regime (E = 0.001). 

We first examine a system where the electron density is small. We take N e = 150 (^ 
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Fig. 10. The LHS and the RHS of eq. (21) plotted against the angular frequency w. The circles 
represent the LHS and the squares represent the RHS. The error bars indicate the standard 
deviations over 20 samples for the LHS, and over 10 samples for the RHS. 



iV a = 4100), iVi = 0, G = 0, K' = 0, G' = 0. In this case many-body effects are expected to be 
small. By measuring the current during 1000 < t < 11000, we get a(T). We investigate a(T) 
at T = 0.125, 0.25, 0.5, 1, for which the lattice doesn't melt according to eq. (7). 
The results are shown in Fig. 11. The solid line is a fitting curve 

a{T) = aT b (22) 

with a = 3.27(±0.38), b = -1.62(±0.09). The observed exponent -1.62 is explained roughly 
as follows. If we assume that the long wave length approximation is valid, the scattering cross 
section satisfies s qU antum(r) oc T. Since s classical of our model imitates s quant um, 

l(T) 1 



< T ) « r(T) * c — , (23) 

holds in the linear response regime, where r is the mean free time, and I is the mean free path, 



yj {v 2 ) T is the standard deviation of the speed of the electrons at temperature T, and n a is 
the density of atoms. n a is independent of T, since the thermal expansion is negligible in this 
model. Therefore, a(T) oc T _1 T C ~ T~ 3 / 2 , which is almost consistent with our result. 

We also analyze the case where the electron density is moderate, N c = 1500 (~ iV a = 
4100). We measure the current during 2000 < t < 3000. The results are shown also in Fig. 11. 
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The dotted line is a fitting curve of eq. (22) with a = 4.31(±0.20), b = -1.41(±0.04). This 
exponent b = — 1.41 ± 0.04 is slightly different from the above case (6 = —1.62 ± 0.09). This 
difference arises from electron-electron interactions, which are treated on an equal footing 
with electron-phonon interactions in our model. 

Those results for the exponent b are consistent with experimental results within experimen- 
tal errors. 22 ~ 24 ) We have thus obtained reasonable dependence of the electrical conductivity 
on temperature, although our model is a classical model. This suggests that we would be able 
to get the correct temperature dependence of the electrical conductivity also in the nonlinear 
nonequilibrium regime. 




temperature T 

Fig. 11. Temperature dependence of the electrical conductivity for N e = 150 (squares) and for N c = 
1500 (circles). The error bars indicate the standard deviations over five impurity distributions. 
The solid, and dotted lines are the fitting curves, eq. (22), with a = 3.27, b = —1.61 for N c = 150, 
and a = 4.31, b = —1.41 for iV = 1500, respectively. 

6. Summary and Conclusion 

In summary, we have proposed a model of electrical conductors, which models interacting 
many electrons scattered by lattice vibrations and random potentials. The model is suitable 
for an MD simulation. 

The model consists of classical electrons and atoms. The atoms compose a lattice vibration 
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system. The electrons are scattered by electron-electron, electron-atom and electron-impurity 
interactions. We have argued that the scattering cross section is physically more important 
than the functional form of a scattering potential. Accordingly, we have set the electron-atom 
interaction potential in such a way that its scattering cross section agrees with that of quantum 
mechanical one. 

To illustrate the advantages of the proposed model, we have applied the proposed model to 
a doped semiconductor at room or higher temperature, and have performed an MD simulation. 

In the linear nonequilibrium regime, we have confirmed the dispersion relations and the 
fluctuation-dissipation relation. We have also obtained reasonable dependence of the electrical 
conductivity on temperature, despite the fact that our model is a classical model. 

In the nonlinear nonequilibrium regime, NESSs are well realized, in which the current I 
and the kinetic energy K fluctuate around the average values with relatively small magnitudes. 
For large E, the response of / to E becomes strongly nonlinear, and the electron temperature 
raises nonlinear ly as a function of the Joule heat. 

Because of these realistic physical properties, we expect that the present model should be 
a good model to explore nonequilibrium states of electrical conductors by an MD simulation. 

Although we have studied macroscopically homogeneous conductors in this paper, we can 
also study macroscopically inhomogeneous conductors using the present model. For example, 
we can analyze conductors that are driven by reservoirs, where the electrical current is induced 
by the chemical potential difference of the reservoirs attached to the both ends. Such systems 
are suitable for exploring nonequilibrium states, because nontrivial and important effects such 
as nonmechanical forces appear, 11,12 ) and may be used to confirm the scaling law of the excess 
noise. 25 ) 
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